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Abstract 

Background: Phytophthora root and stem rot (PRR) of soybean, caused by Phytophthora sojae, is controlled by Rps 
genes. However, little is known regarding the /?p5-induced molecular responses to P. sojae and how they actually 
overlap. We thus sequenced, analyzed, and compared the transcriptomes of 10 near isogenic lines (NILs), each with 
a unique Hps gene/allele, and the susceptible parent Williams, pre- and post-inoculation with the pathogen. 

Results: A total of 4,330 differentially expressed genes (DEGs) were identified in Williams versus 2,014 to 5,499 DEGs 
in individual NILs upon inoculation with the pathogen. Comparisons of the DEGs between the NILs and Williams 
identified incompatible interaction genes (HGs) and compatible interaction genes (CIGs). Hierarchical cluster and 
heatmap analyses consistently grouped the NILs into three clusters: Cluster I {Rpsl-a), Cluster II {Rpsl-b, 1-c and 1-k) 
and Cluster III (Rps3-a, 3-b, 3-c, 4, 5, and 6), suggesting an overlap in /?p5-induced defense signaling among certain 
NILs. Gene ontology (GO) analysis revealed associations between members of the WRKY family and incompatible 
reactions and between a number of phytohormone signaling pathways and incompatible/compatible interactions. 
These associations appear to be distinguished according to the NIL clusters. 

Conclusions: This study characterized genes and multiple branches of putative regulatory networks associated with 
resistance to P. sojae in ten soybean NILs, and depicted functional "fingerprints" of individual /?p5-mediated 
resistance responses through comparative transcriptomic analysis. Of particular interest are dramatic variations of 
detected DEGs, putatively involved in ethylene (ET)-, jasmonic acid (JA)-, (reactive oxygen species) ROS-, and 
(MAP-kinase) MAPK- signaling, among these soybean NILs, implicating their important roles of these signaling in 
differentiating molecular defense responses. We hypothesize that different timing and robustness in defense 
signaling to the same pathogen may be largely responsible for such variations. 
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Background 

Phytophthora root and stem rot (PRR) is one of the 
most devastating diseases of soybean (Glycine max), 
causing nearly $200 million in annual yield losses in the 
U.S. alone [1]. PRR is caused by the soil-borne, hemibio- 
trophic pathogen Phytophthora sojae, and is most effect- 
ively controlled by Rps (Resistance to P. sojae) genes. 
The resistance conferred by an Rps gene is race specific, 
and the interaction between an Rps gene and the corre- 
sponding avirulence (Avr) gene in P. sojae, follows the 
gene-for-gene model [2]. Currently, 18 Rps genes/alleles 
from soybean and 12 Avr genes from P. sojae have been 
identified [3-5]. 

Like most resistance (R) genes, the Rps gene family en- 
codes, or is predicted to encode nucleotide binding leu- 
cine rich repeat (NB-LRR)-type proteins [6-8], which are 
able to recognize the Avr effector proteins of pathogens 
directly or indirectly to induce the appropriate defense 
response [9,10]. The first evidence of a direct interaction 
between an R and Avr protein was reported in the flax- 
Melampsora lini pathosystem [11]. Indirect interactions, 
however, are more common. In these cases, the R pro- 
tein appears to require the existence of a 'guard protein' 
or 'decoy' in order to recognize an Avr protein [12,13]. 
A classical example of this type of interaction is the 
RPM1 -interacting protein 4 (RIN4), which is 'guarded' 
by RPM1 (RESISTANCE TO PSEUDOMONAS SYRIN- 
GAE PV MACULICOLA 1) and RPS2 (RESISTANCE 
TO P. SYRINGAE 2) proteins in Arabidopsis [14,15]. 

The recognition of pathogen effectors triggers an in- 
nate immunity response that is mediated by two distinct 
types of NB-LRR proteins, the toll-interleukin-1 receptor 
(TIR)-NB-LRR proteins and the coiled-coil (CC)-NB- 
LRR proteins. The former requires EDS1 (ENHANCED 
DISEASE SUSCEPTIBILITY 1), a central regulator of 
effector-triggered immunity (ETI), which functions to- 
gether with PAD4 (PHYTOALEXIN DEFICIENT 4). In 
contrast, the latter NB-LRR proteins are independent of 
EDS1 but require NDR1 (NON-RACE-SPECIFIC DIS- 
EASE RESISTANCE 1) [16-21]. The interaction among 
these intracellular proteins results in a regulation network 
of phytohormone signaling, which is mainly mediated by 
salicylic acid (SA) for biotrophic and hemibiotrophic path- 
ogens, and jasmonic acid (JA) and ethylene (ET) for 
necrotrophic pathogens [8,22]. In addition to SA, JA and 
ET, other phytohormones such as brassinosteroids (BR), 
abscisic acid (ABA), auxin, gibberellins (GA) and cytokinin 
(CK) also contribute to plant immune responses, with 
complex crosstalk between one another [23]. 

It is suggested that the resistance to P. sojae conferred 
by Rps genes is mediated by the SA signaling pathway, 
with the induction of SA-mediated systemic acquired re- 
sistance (SAR) occurring several days post infection 
(dpi) via expression of the gene NPR1 [24-26]. During 



this process, two putative regulators of the chromosome 
condensation 1 (RCC1) gene family are down-regulated dur- 
ing the incompatible interaction [27]. Besides SA, exogenous 
treatment of 1-aminocyclopropane-l-carboxylic acid (ACC, 
a precursor of ET) has been shown to enhance resistance 
while applications of GA and ABA induce susceptibility, 
highlighting the complexity of phytohormone signaling 
pathways in response to attack by P. sojae [26,28,29]. 

A microarray study of transcriptomes from one sus- 
ceptible and two partially resistant soybean genotypes 
indicated that 97-99% of all detectable genes experienced 
transcriptional modulation five dpi in response to infec- 
tion by P. sojae [30]. However, the majority of these dif- 
ferences were less than two fold. Another microarray 
study of transcriptomic changes in soybean revealed a 
peak in most defense related genes at 24 hours post in- 
oculation (hpi) with P. sojae [24]. Recently, Zhang et al 
conducted a proteomics study in which 46 differentially 
expressed proteins were identified in soybean after infec- 
tion with P. sojae [31]. Among these, 26 were affected 
during the incompatible interaction, while the other 20 
were altered during the compatible interaction. These 
studies have contributed to our understanding of the 
interaction between soybean and P. sojae. Nevertheless, 
what mechanisms underlie compatible and incompatible 
interactions, how molecular responses mediated by a 
variety of Rps genes/alleles resemble or differ from one 
another, as well as the nature of these similarity or dif- 
ference, remain largely unknown. 

Access to the soybean genome sequence of Williams 82 
(Rpsl-k) [32], and the advent of high-throughput RNA se- 
quencing (RNA-Seq) by next-generation sequencing tech- 
nology, have allowed researchers to take a novel look at the 
molecular interaction between soybean and P. sojae. Kim 
et al reported the first application of RNA-Seq for profiling 
gene expression in soybean in response to pathogen attack 
[33]. In this study, the transcriptomes of two near isogenic 
lines (NILs), one resistant and one susceptible to bacterial 
leaf pustule, were analyzed 0, 6, and 12 hpi and a total of 
2,761 differentially expressed genes (DEGs), including a set 
of defense response genes, were identified. 

NILs are pure breeding lines that are developed by 
backcrossing a donor line with the recurrent parent for 
at least five generations to achieve introgression of a de- 
sired trait. As such, they share more than 98% genetic 
identity with the recurrent parent, except for the region 
where the desired gene is located. In order to determine 
the molecular mechanisms responsible for mediated 
defense and to understand molecular basis for diverse 
responses to the pathogen, we analyzed the transcrip- 
tomes of 10 NILs, each with a unique Rps gene/allele, 
along with the susceptible recurrent parent Williams, 
pre- and post-inoculation with a race 1 isolate of P. sojae 
(avirulent towards NILs; virulent towards Williams). 
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Results 

Symptom development in 10 NILs and the recurrent 
susceptible parent inoculated with P. sojae 

The 10 soybean NILs carrying Rpsl-a, Rpsl-b, Rpsl-c, 
Rpsl-k, Rps3-a, Rps3-b, Rps3-c, Rps4, RpsS, or Rps6, 
within the Williams background, provided a unique re- 
source for investigating the common and specific defense 
responses mediated by individual Rps genes against P. 
sojae (Additional file 1). Phenotypic reactions of the 10 
NILs and Williams to P. sojae 7 days post-inoculation 
(dpi) with the pathogen were evaluated under greenhouse 
conditions. In each of three experimental replications, ap- 
proximately half the number of seedlings from each line 
was inoculated with race 1 of P. sojae. The remaining 
seedlings were "mock" inoculated in the same manner 
without the pathogen. In each experimental replicate, Wil- 
liams was susceptible (expanding lesion/plant death) to P. 
sojae, while all NILs were resistant (Additional file 2). NILs 
containing Rpsl-a, Rpsl-b, Rpsl-c, Rpsl-k and Rps4 dis- 
played 100% survival when challenged with the pathogen. 
In contrast, NILs containing Rps3-a, Rps3-b, Rps3-c, RpsS 
and Rps6 showed a slight variation in the proportion of 
surviving seedlings across the three replicates, which may 
be attributed to minor differences in environmental condi- 
tions. Despite this variation, the percentage of surviving 
seedlings in each replicate was equal to or greater than 
75%, which is generally used as a criterion for defining a 
pure line resistance [34]. All mock-inoculated seedlings 
were asymptomatic of PRR (Additional file 2). 

Transcriptional changes in 10 NILs and the recurrent 
susceptible parent in response to P. sojae 

RNAs, representing the entire transcrip tomes of P. sojae 
inoculated and mock-inoculated seedlings from each of 
the 10 NILs and Williams 24 hpi with the pathogen, 
were sequenced and 14.5 to 50.1 million raw reads were 
generated for individual samples (Table 1). Trimming 
adaptor sequences and elimination of low quality reads 
and reads shorter than 101 bp, resulted in 13.5 to 39.5 
million clean reads for each sample. Among these, 71.5% 
to 87.8% per sample were uniquely mapped to the soy- 
bean reference genome vl.O (Table 1). 

Based on the 46,367 high-confidence gene models anno- 
tated in the soybean reference genome [32], the relative 
change in the abundance of reads mapped to the same 
genes for inoculated and mock-treated samples, and the 
criteria for defining DEGs pre- and post-inoculation, a 
total of 9,847 non-redundant gene models were identified 
as DEGs in at least one of the 11 soybean lines. The num- 
bers of DEGs identified between Williams and the 10 NILs 
or among the 10 NILs vary greatly. For example, 2,014 to 
2,366 DEGs were found in individual NILs containing 
Rpsl-a, Rpsl-b, Rpsl-c, or Rpsl-k, while 3,038 to 5,499 
DEGs were detected in individual NILs containing Rps3-a, 



Rps3-b, Rps3-c, Rps4, RpsS or Rps6 (Table 2). Among the 

10 NILs, 1,274 to 2,823 DEGs were up-regulated, while 
643 to 2,744 DEGs were down-regulated in response to in- 
oculation with the pathogen (Table 2). In comparison, 
2,460 up-regulated and 1,870 down- regulated DEGs were 
identified in Williams upon inoculation with P. sojae. 

To validate DEGs profiled by RNA-Seq analysis, six de- 
tected DEGs per soybean line (Additional file 3) were ran- 
domly chosen for quantitative RT-PCR (qRT-PCR) analysis 
with the same RNA samples as used for RNA-Seq. Sub- 
sequently, the patterns of differential expression of these 
genes detected by RNA-Seq and qRT-PCR were compared. 
As shown in Additional file 4, significant correlations be- 
tween the patterns of DEGs detected by the two approaches 
were observed for each set of genes chosen from 11 indi- 
vidual lines (Pearsons correlation coefficient (r), ranging 
from 0.725 to 0.994), indicating the reliability of DEGs 
identified by RNA-Seq (Additional file 4). 

Characterization of genes associated with incompatible, 
compatible, and basal interactions 

In an attempt to decipher the molecular basis of resist- 
ance and susceptibility to P. sojae, we applied a com- 
parative transcrip tomics approach to identify DEGs 
specifically associated with a host response. Three major 
groups of genes were classified: 1) Incompatible inter- 
action genes (IIGs), are DEGs identified in NILs and as- 
sociated with resistance; 2) Compatible interaction genes 
(CIGs) are DEGs identified in Williams and specifically 
associated with the susceptible response; and 3) Basal 
interaction genes (BIGs), which are DEGs shared by all 
NILs and Williams (Figure 1). 

A total of 5,806 non-redundant IIGs, 1139 CIGs, and 
835 BIGs were identified (Figure 1). The number of up- 
regulated and down-regulated IIGs range from 159 to 
1,158 and from 141 to 2,017, respectively, among the 10 
NILs. Of the 1,139 CIGs, 369 were up-regulated and 770 
were down-regulated; and of the 835 BIGs, 696 were up- 
regulated and 139 were down-regulated (Figure 1). 

Clusters of /?ps-mediated IIGs 

To understand how transcriptomic changes mediated by 
different Rps genes/alleles overlap, we performed hierarch- 
ical cluster analysis using pvclust, an R package for asses- 
sing the uncertainty in hierarchical clustering [35], with 
log 2 fold change (FC) of the 5,806 IIGs identified in at 
least one of the ten NILs. The NILs were grouped into 
three clusters, Clusters I, II, and III (Figure 2). Cluster I 
(C-I) consisted of the NIL containing Rpsl-a only. Cluster 

11 (C-II) was composed of those NILs containing Rpsl-b, 
Rpsl-c and Rpsl-k, while Cluster III (C-III) included those 
NILs containing Rps3-a, Rps3-b, Rps3-c, Rps4, RpsS, and 
Rps6. This clustering of IIGs was further supported by the 
hierarchical cluster structure of the ten NILs revealed by 
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Table 1 Statistics of pair-end reads in RNA sequencing experiment 
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a Reads with a quality score less than 20 and length less than 101 nucleotides were excluded. 
b Reads uniquely mapped to G. max genome assembly (vl.0) using TopHat 2.0 software [58]. 
Percentage of uniquely mapped reads out of clean reads. 



heatmap analysis [36] (Additional file 5). In addition, the 
numbers of IIGs identified for each individual NIL, were 
unevenly distributed among different clusters (Figures 1 
and 2). For example, 300, 1,193-3,073, and 556-614 IIGs 
were found in C-I, C-II, and C-III, respectively. The 

Table 2 Differentially expressed genes (DEGs) in each 
soybean line 



Soybean line 


Up-regulated 
DEGs 


Down-regulated 
DEGs 


Total 


Williams ifps) 


2460 


1870 


4330 


Union (Rpsl-a) 


1432 


643 


2075 
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4613 
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2755 


2744 


5499 
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2033 


1005 


3038 



numbers of IIGs in NILs containing different alleles at a 
same gene locus also varied (Figures 1 and 2). 

The IIGs among the three clusters were analyzed by both 
broad range comparison, where all the genes within a clus- 
ter were counted, and narrow range comparison, where 
only those genes shared by the NILs within a cluster were 
counted. Broad range comparison showed that 75 up- 
regulated and 68 down-regulated IIGs were shared by the 
three clusters, whereas 62, 214, and 4,488 IIGs were unique 
to NILs in Clusters I, II, and III, respectively (Figure 3A,B). 
In contrast, narrow range comparison indicated that only 
16 up-regulated and 7 down-regulated IIGs were shared 
among the three clusters (Figure 3C,D). The annotation of 
these 23 IIGs is listed in Table 3. 



Putative functions of DEGs based on gene ontology (GO) 
analysis 

To understand the functional components involved in 
molecular responses to P. sojae, we annotated the puta- 
tive products encoded by IIGs, CIGs, and BIGs based on 
GO analysis (Figure 4; Additional file 6). These DEGs 
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Figure 1 Comparison of differentially expressed genes (DEGs) between Williams (susceptible to Phytophthora sojae) and 10 near 
isogenic lines (NILs) each containing a single Rps (Resistance to P. sojae) gene. A: up-regulated DEGs. B: down-regulated DEGs. Red and 
purple represent the number of DEGs specific to individual Rps genes and collectively referred to as incompatible interaction genes. Brown and 
dark purple represent the number of DEGs shared by an individual NIL and Williams. Green and light blue represent the number of DEGs 
specifically expressed in 'Williams' when compared to a specific NIL. The central green and light blue circles represent the common proportion 
of Williams specific DEGs that are collectively referred to as compatible interaction genes. 
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Figure 2 Hierarchical cluster analysis (pvclust) of incompatible 
interaction genes identified for 10 soybean near isogenic lines, 
each containing a single Rps gene. The numbers behind each 
node represent AU (Approximately Unbiased) / BP (Bootstrap 
Probability) for an estimation of the confidence of each node. 



were grouped into seven categories: 1) Response to bi- 
otic or abiotic stress; 2) Biological regulation; 3) Growth 
and development; 4) Transport; 5) Metabolism; 6) Mis- 
cellaneous; and 7) Unclassified. For IIGs, the majority of 
annotated genes fell into the category of "response to bi- 
otic stress and abiotc stress" (Figure 4; Additional files 6 
and 7). This pattern was also observed for CIGs and 



B 



12 



C-l 

35 

75 



37 



C-l 
27 

68 



43 



C-ll 

79 



385 



C-lll 
1481 



C-ll 
135 



471 



C-lll 

3007 



C-l 
123 

16 



C-l 
120 



13 



11 



C-ll 
13 



44 



C-lll 

282 



C-ll 
16 



15 



C-lll 
158 



Figure 3 Comparison of incompatible interaction genes 
between three clusters, C-l, C-ll, and C-lll. A: Broad range 
comparison for up-regulated genes. B: Broad range comparison for 
down-regulated genes. C: Narrow range comparison for up- 
regulated genes. D: Narrow range comparison for down- 
regulated genes. 
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Table 3 Incompatible interaction genes (HGs) shared by all the near isogenic lines and their annotation 



Gene 


Arabidopsis 
homologous gene 


Gene symbol 


Gene annotation 


GO group 


Up regulated 










Glyma01g32130 


AT2G40435 


- 


Transcription factor, bHLH 


Response to biotic or abiotic stress 


Glyma02g16710 


AT1G03220 


- 


Eukaryotic aspartyl protease family protein 


Metabolism 


Glyma03g01010 


AT4G26010 


- 


Peroxidase superfamily protein 


Response to biotic or abiotic stress 


Glyma05g34870 


AT1G 14870 


PCR2 


Plant cadmium resistance 2 


- 


Glyma06g09220 


AT5G25880 


NADP-ME3 


NADP-malic enzyme 3 


Response to biotic or abiotic stress 


Glyma06g 13090 


AT5G24110 


WRKY30 


WRKY DNA-binding protein 30 


Response to biotic or abiotic stress 


Glyma09g 15090 


AT4G21380 


RK3 


Receptor kinase 3 


Metabolism 


Glyma11g07430 


AT5G03260 


LAC11 


Laccase 1 1 


Metabolism 


Glyma13g17800 


AT4G04450 


WRKY42 


WRKY family transcription factor 


Biological regulation 


Glyma13g30770 


AT1G28480 


GRX480 


Thioredoxin superfamily protein 


Response to biotic or abiotic stress 


Glyma13g36070 


AT3G20660 


4-Oct 


Organic cation/carnitine transported 


Biological regulation 


Glyma14g39300 


AT1G66160 


CMPG1 


CYS, MET, PRO, and GLY protein 1 


Response to biotic or abiotic stress 


Glyma15g05810 


AT2G41480 


- 


Peroxidase superfamily protein 


Response to biotic or abiotic stress 


Glyma15g06010 


AT5G64260 


EXL2 


Exordium like 2 


- 


Glyma15g 19600 


AT2G05940 


- 


Protein kinase superfamily protein 


Metabolism 


Glyma20g38000 


AT1G09090 


RBOHB 


Respiratory burst oxidase homolog B 


Response to biotic or abiotic stress 


Down regulated 










Glyma04g 12290 


AT1G71380 


CEL3 


Cellulase 3 


Response to biotic or abiotic stress 


Glyma09g37290 


AT5G 15230 


GASA4 


GAST1 protein homolog 4 


Response to biotic or abiotic stress 


Glyma12g 10670 


AT3G21700 


SGP2 


Ras-related small GTP-binding family protein 


Response to biotic or abiotic stress 


Glyma12g29980 


AT4G39370 




Ubiquitin-specific protease 27 




Glyma13g44210 


AT5G20190 




Tetratricopeptide repeat (TPR)-like superfamily protein 




Glyma15g 18440 


AT4G29720 


PA05 


Polyamine oxidase 5 


Metabolism 


Glyma18g49400 


AT5G 15230 


GASA4 


GAST1 protein homolog 4 


Response to biotic or abiotic stress 



BIGs (Figure 4; Additional file 7). Of the 23 IIGs shared 
by the 10 NILs, 12, 2, and 5, were found to be related to 
"response to biotic or abiotic stress", "biological regula- 
tion", and "metabolism", respectively (Table 3). The 
remaining four shared genes could not be classified into 
any functional category. 

Putative transcription factors (TFs) involved in molecular 
responses to P. sojae 

It is documented that TFs play important roles in plant 
defense and stress responses. As such, we annotated and 
analyzed putative TFs in the sets of DEGs detected in this 
study based on a soybean TF database that is composed of 
5,671 predicted TFs [37]. Of the 5,806 non-redundant IIGs, 
282 up-regulated and 543 down-regulated genes were an- 
notated as putative TFs. The up-regulated and down- 
regulated IIG TFs range from 13 to 177 and from 21 to 
307, respectively, among the 10 NILs (Additional files 8 
and 9). The number of IIG TFs also varied among the three 
clusters, with 34, 164, and 791 in C-I, C-II, and C-III, 



respectively. The most abundant ones along these IIG TFs 
(either up-regulated or down-regulated) were of the WRKY 
family, accounting for 23.1% of all IIG TFs identified in C-I, 
18.5-25.8% in C-II, and 14.0-18.6% in C-III (Additional file 
9). Only three up-regulated IIG TFs (1 bHLH and 2 
WRKY) and one down-regulated TF (TPR) were found to 
be shared by all the 10 NILs (Table 3; Additional file 9). 

Of the 1,139 CIGs identified in Williams, 43 up- 
regulated and 149 down-regulated genes were found to 
be putative TFs (Additional file 9). The largest propor- 
tion of up-regulated CIG TFs belonged to the MYB/HD- 
like family (16.3%), followed by the AP2-EREBP family 
(14.0%). In contrast, members of the AP2-EREBP family 
represented the largest number (20.8%) of down-regulated 
CIG TFs, followed by MYB/HD-like (10.1%) and bHLH 
(9.4%). None of these 43 up-regulated CIG TFs belonged 
to WRKY family. 

Of the 835 non-redundant BIGs shared by all the 11 
soybean lines, 41 up-regulated and 9 down- regulated 
ones were found to be TFs (Additional file 9). Among 
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Figure 4 Heatmap analysis of three categories of differentially expressed genes based on gene ontology analysis. The values used to 
draw heatmaps are Log 2 (fold change) of expression level of post inoculation to mock inoculation. A: Incompatible interaction genes. The Log 2 (fold 
change) less than 3 in all 10 resistant near isogenic lines were not shown here. B: Compatible interaction genes. C: Core basal interaction genes. 



the up-regulated BIG TFs, the most abundant ones are 
of the WRKY family (34%), followed by the MYB/HD- 
like family (19.5%). 

Knowledge-based comparative analysis of ftps 
gene-mediated defense signaling pathways 

The DEGs showing homologies to previously identified 
protein genes responsible for pathogen recognition and 
defense, or defense- related phytohormone signaling genes 
were further analyzed. Overall, these DEGs exhibited di- 
verse patterns of distribution among the three NIL clus- 
ters (Figure 5). For example, of the 26 up-regulated IIGs 
homologous to previously reported defense/stress signal- 
ing genes, 24 were found only in NILs within Cluster III. 
In contrast, only one (Glyma20g38000) of the two RBOH 
B gene homologs was found to be up-regulated in all 10 
NILs. It is notable that the majority of the putative DEGs 
involved in defense-related phytohormone signaling path- 
ways showed distinct or opposite patterns of variation in 
gene expression between Williams and NILs, particularly 
the NILs within the Cluster III (Figure 5). 

Previous studies indicated that the resistance to P. 
sojae in soybean was mediated by SA signaling [24-26], 



with the NPR1 as the key component of SA-mediated 
signaling [38]. In this study, two NPR1 -like IIGs 
(Glyma02g4S260 and Glymal4g03510) were identified 
(Figure 5). However, these two genes differ from the two 
NPR1 -like genes reported by Sandhu et al [25]. In the later 
study, they studied soybean leaves four dpi with P. sojae 
race 4 while, in the present study, we studied stems 24 hpi 
with P. sojae race 1. As such, these two studies are less 
comparable. 

As an antagonist of SA, the JA pathway is repressed by 
JAZ proteins in Arabidopsis [39,40]. Several JAZ homologs 
were identified as IIG and/or CIG, including homologs of 
JAZ1, JAZ6, JAZ8 and JAZ12 (Figure 5). These homologs 
showed a consistent expression pattern that was up- 
regulated NILs and down-regulated in Williams. In contrast 
to the JAZ proteins, JAR1 (Jasmonate resistantl) appears to 
be required for the activation of the JA signaling pathway 
in Arabidopsis [41,42]. We found that a soybean JAR1 
homolog (Glyma07g06370) was up-regulated in Williams 
upon inoculation with the P. sojae. These observations sug- 
gest that the JA signaling pathway may have played oppos- 
ite roles in phytohormone signaling between incompatible 
and compatible interactions. 
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GeneJD Gene 
Symbol 


Wms 


Rps1-a 

C-I b 


Rps1-b 

C-ll 


Rps1-c 

C-ll 


Rps1-k 

C-ll 


Rps3-a 
C-lll 


Rps3-b 
C-lll 


Rps3-c 
C-lll 


Rps4 
C-lll 


Rps5 
C-lll 


Rps6 
C-lll 


Gene 
Group 3 


Glyma04g38700 PAD4 


-0.79 


-0.59 


-0.69 


-0.66 


-0.78 


-1.29 


-0.93 


-1.08 


-1.32 


-1.25 


-1.08 


IIG 


Glyma05g32760 PAD4 


-1 .27 


-0.48 


0.12 


0.04 


0.53 


0.91 


0.91 


0.68 


1.04 


1.13 


0.91 


CIG/IIG 


Glyma08g00420 PAD4 


-0.47 


0.11 


0.23 


0.19 


0.66 


1.21 


1.17 


1.16 


1.17 


1.41 


1.08 


IIG 


Glyma12g34210 NDR1 


-0.49 


-0.43 


0.39 


0.62 


0.06 


0.84 


1.45 


0.87 


0.76 


1.02 


1.05 


IIG 


Glyma03g 19920 RIN4 


0.74 


0.67 


0.46 


0.59 


-0.01 


1.21 


0.99 


1.50 


1.22 


1.42 


0.94 


IIG 


Glyma08g08810 FLS2 


-0.95 


-0.11 


-0.45 


-0.56 


0.48 


0.86 


1.23 


2.09 


1.16 


1.23 


1.66 


IIG 


Glyma13g20800 PAL 


0.90 


0.51 


0.61 


0.63 


0.46 


0.76 


0.81 


1.08 


1.01 


1.19 


1.01 


IIG 


Glyma03g39610 RBOH B 


-0.91 


-0.15 


0.85 


1.23 


1.22 


2.52 


3.13 


1.73 


2.25 


2.69 


1.92 


IIG 


Glyma20g38000 RBOH B 


0.93 


1.40 


1.48 


1.95 


1.74 


1.83 


2.11 


2.27 


1.78 


1.65 


1.42 


IIG 


Glyma11g15700 MPK3 


-0.52 


-0.18 


0.26 


0.22 


0.29 


1.33 


1.75 


1.14 


1.44 


1.64 


1.10 


IIG 


Glyma12g07770 MPK3 


-0.40 


-0.10 


0.36 


0.54 


0.44 


1.50 


1.99 


1.45 


1.37 


1.86 


1.36 


IIG 


Glyma07g07270 MPK4 


0.68 


0.62 


1.45 


1.31 


0.45 


2.01 


2.11 


2.15 


2.23 


2.39 


2.09 


IIG 


Glyma09g39190 MPK4 


0.93 


0.43 


0.59 


0.10 


0.17 


0.59 


0.86 


1.48 


0.93 


1.19 


1.20 


IIG 


Glyma08g02060 MPK4 


-0.46 


-0.80 


-0.03 


-1.78 


-5.77 


-1.25 


-0.39 


-0.83 


-1.74 


-1.80 


-0.34 


IIG 


Glyma02g45260 NPRI-like 


0.48 


0.67 


0.34 


0.58 


0.46 


0.99 


1.04 


0.73 


1.23 


0.90 


0.72 


IIG 


Glyma14g03510 NPRI-like 


1.18 


0.75 


1.60 


0.71 


1.08 


0.52 


1.66 


1.10 


-0.02 


1.73 


0.48 


IIG 


Glyma09g08290 JAZ1 


-0.84 


-0.48 


0.23 


0.63 


0.14 


0.95 


1.21 


0.47 


0.82 


1.22 


0.56 


IIG 


Glyma11g04130 JAZ1 


-2.06 


-0.67 


0.67 


0.56 


0.63 


2.15 


2.21 


1.49 


1.88 


2.01 


1.10 


IIG/CIG 


Glyma13g17180 JAZ1 


-1.09 


-0.40 


0.45 


0.48 


0.71 


1.61 


1.93 


1.27 


1.50 


1.84 


1.17 


IIG/CIG 


Glyma15g 19840 JAZ1 


-0.98 


-0.37 


0.33 


0.28 


0.28 


1.00 


1.09 


0.39 


0.78 


1.15 


0.53 


IIG 


Glyma17g05540 JAZ1 


-2.02 


-0.93 


0.24 


0.45 


0.45 


1.59 


1.46 


0.66 


1.54 


1.42 


0.57 


IIG/CIG 


Glyma07g04630 JAZ6 


-0.83 


-0.73 


0.14 


0.20 


0.41 


1.02 


0.92 


0.82 


1.08 


1.24 


0.57 


IIG 


Glyma16g01220 JAZ6 


-0.91 


-0.43 


0.41 


0.42 


0.65 


1.45 


1.37 


1.05 


1.40 


1.32 


0.84 


IIG 


Glyma05g27280 JAZ8 


-1.59 


-1 .72 


2.40 


0.95 


0.32 


3.85 


2.08 


1.68 


3.05 


3.16 


2.73 


IIG 


Glyma08g 10220 JAZ8 


0.09 


-0.55 


1.22 


0.94 


0.97 


3.35 


3.01 


3.12 


3.43 


4.12 


2.94 


IIG 


Glyma13g29070 JAZ8 


-3.00 


-1.08 


0.21 


0.52 


0.50 


2.13 


2.27 


1.70 


2.57 


1.94 


1.28 


IIG/CIG 


Glyma13g 17640 JAZ12 


-0.82 


-0.19 


0.63 


0.27 


0.61 


1.65 


0.35 


0.86 


1.50 


1.03 


0.21 


IIG 


Glyma17g04850 JAZ12 


-1.05 


-0.59 


0.05 


-0.49 


0.42 


1.04 


-0.49 


0.44 


0.07 


0.27 


-0.23 


CIG 


Glyma07g06370 JAR1 


1.11 


0.46 


0.21 


0.68 


0.53 


0.42 


0.10 


0.47 


-0.06 


0.22 


0.39 


CIG 


Glyma20g21780 ETR2 


-0.56 


-0.98 


-0.41 


-0.27 


-1.06 


-0.69 


-0.71 


-0.72 


-0.84 


-1.00 


-0.27 


IIG 


Glyma20g34420 ETR2 


1.33 


0.45 


0.89 


0.84 


0.34 


0.08 


0.41 


0.73 


-0.34 


-0.29 


0.40 


CIG 


Glyma03g41220 EIN4 


-0.47 


-0.68 


-0.90 


-1.13 


-1.37 


-2.65 


-1.86 


-1.88 


-2.05 


-2.40 


-1.04 


IIG 


Glyma19g43840 EIN4 


-0.41 


-0.19 


-0.78 


0.00 


-0.40 


-0.94 


-0.60 


-0.33 


-1.38 


-1.44 


-0.18 


IIG 


Glyma04g07110 EBF1 


0.79 


-0.08 


-0.16 


0.62 


-0.91 


-0.77 


-0.68 


-0.94 


-1.39 


-0.79 


-0.56 


IIG 


Glyma04g20330 EBF1 


2.68 


0.79 


0.29 


1.04 


-0.26 


0.16 


0.52 


0.35 


-0.51 


-0.08 


0.98 


CIG 


Glyma06g07200 EBF1 


1.18 


-0.02 


-0.27 


0.88 


-0.79 


-0.53 


-0.58 


-0.75 


-1.29 


-0.88 


-0.22 


IIG/CIG 


Glyma13g23510 EBF1 


1.71 


0.37 


0.38 


0.94 


-0.32 


0.04 


0.13 


-0.49 


-0.79 


-0.38 


0.57 


CIG 


Glyma17g31940 EBF1 


1.32 


0.01 


-0.24 


0.95 


-0.98 


-0.95 


-0.83 


-0.83 


-1.56 


-0.99 


-0.36 


IIG/CIG 


Glyma05g24770 BAK1 


0.19 


1.14 


2.02 


2.40 


2.03 


2.88 


3.19 


2.14 


2.84 


3.40 


2.27 


IIG 


Glyma08g 19270 BAK1 


0.75 


0.67 


0.94 


0.76 


0.68 


0.75 


1.22 


1.16 


0.92 


0.85 


0.76 


IIG 


Glyma15g05730 BAK1 


0.85 


0.82 


0.75 


0.84 


0.77 


0.95 


1.07 


1.43 


0.99 


0.98 


0.93 


IIG 



Figure 5 Expression of selected MGs or CIGs involved in defense response and defense related phytohormone signaling pathways. The 

values and the corresponding depth of different colors in each square indicate Log 2 (fold change) of expression level of post inoculation to mock 
inoculation. Classification of differentially expressed genes. IIG = Incompatible interaction gene. CIG = Compatible interaction gene. Clustering of 
Rps genes based on hierarchical cluster analysis of IIGs. C-l = Cluster I. C-ll = Cluster II. C-lll = Cluster III. 



The ET signaling pathway is generally considered to 
work cooperatively with JA in response to necrotrophic 
pathogens [8,22], In Arabidopsis, ET receptors ETR2 
[43] and EIN4 [44] were both negative regulators of ET 
signaling, and an ubiquitin/proteasome pathway for the 
degradation of EIN3 was mediated by an F-box protein, 
EBF1, to negatively regulate ET responses [23]. We 
found that one homolog of ETR2 and two homologs of 
EIN4 were down- regulated in NILs, while another 
homolog of ETR2 was up-regulated in Williams. In 
addition to ETR2 and EIN4, we identified five EBF1 



homologs that were down-regulated in all NILs and up- 
regulated in Williams (Figure 5), suggesting that the ET 
signaling response to the pathogen was repressed during 
compatible interactions but activated during incompat- 
ible interactions. 

In addition to the three major phytohormones de- 
scribed above, BR have also been found to enhance dis- 
ease resistance against viral and bacterial pathogens in 
tobacco and rice [45]. BAK1 is a crucial component of 
BR signaling, and can interact with LRR-RLK (leucine- 
rich repeat receptor-like protein kinase) genes such as 
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FLS2 to promote their function in plant defense re- 
sponses [46]. As shown in Figure 5, three BAK1 homo- 
logs were identified as IIGs and all were up-regulated 
upon inoculation with the pathogen, suggesting a posi- 
tive role of BR signaling in defense response to P. sojae. 

Discussion 

The nature of distinction in molecular immune response 
among NILs 

One of the most striking observations in this study was 
the distinction in defense response to inoculation with 
P. sojae among the ten resistant NILs as revealed by 
comparative transcriptomic analysis. This was reflected 
not only by the level of variation in the number of IIGs 
(ranging from 300 to 3073) detected in individual NILs, 
but also by the relative small numbers of up-regulated 
(16) and down-regulated (7) IIGs shared by all NILs 
(Table 3). In general, different Rps genes recognize dis- 
tinct Avr genes in a P. sojae population [47], but the Avr 
determinants of different Rps genes in a same isolate of 
the pathogen (i.e., race 1 in this study) may not be dis- 
tinct. One hypothesis that may explain the observed 
variations of DEGs among different NILs may be differ- 
ential timing and robustness in defense responses and 
signaling mediated by different Rps genes/alleles [48]. 
Although the timing of gene expression mediated by dif- 
ferent R genes in response to the same pathogen has not 
been characterized in plants, dramatic and rapid changes 
in gene expression after inoculation with, or following 
infection by a pathogen, have been observed [24,49]. A 
recent study demonstrated that -5% of DEGs were 
shared by soybean lines resistant to bacterial leaf pustule 
at 6 and 12 hpi with Xanthomonas axonopodis [33]. In 
this study, the majority (>95%) of DEGs showed less 
than an eight-fold difference in expression levels upon 
inoculation with the pathogen when given the same cut- 
off time for each of the 11 lines. Such a low level of ex- 
pressional changes and uniformed cutoff may have affected 
the detectability of shared DEGs among different lines. In 
addition, the NILs remain 1-2% genomic difference (mostly 
surrounding the Rps loci) from each other, which may have 
affected expression patterns of a small proportion of genes. 
Furthermore, inoculated regions of seedling stems, near 
which layers of cells respond to the pathogen, vary in size 
among different plants and NILs. Such a variation may 
affect effective identification of DEGs. Therefore, the num- 
bers of DEGs, and particularly IIGs specific to each NIL 
could be over-estimated. Nevertheless, the DEGs in each 
line were detected pre-and post-inoculation with the same 
pathogen, and thus the influence of genomic difference sur- 
rounding Rps genes on our pipeline for detection of DEGs 
triggered by the pathogen may be minimal. 

It is documented that resistance genes recognize path- 
ogens either directly or indirectly by guarding a protein 



or using a decoy [13]. In Arabia 1 opsis, RIN4 is 'guarded' 
by resistance genes RPM1 and RPS2, which trigger the 
immune response [14,15]. In soybean, two of four RIN4- 
homologous genes appear to function as a heteromeric 
complex in mediating RPG1-B and RPM1 -derived resist- 
ance to Pseudomonas syringae, and silencing GmRIN4a 
or GmRIN4b in rpgl-b plants enhances basal resistance 
to virulent strains of P. syringae and the oomycete P. 
sojae [50]. In our study, a single RIN4 homolog was 
found to be up-regulated in C-III NILs, but not in NILs 
within C-I and C-II. If this RIN4 homolog did encode 
the "guard" protein that was recognized by the Rps genes 
in the C-III NILs, then the relatively low numbers of 
IIGs identified among NILs in C-I and C-II, as compared 
with the C-III NILs, may be explained by the distinct 
"recognition" processes needed for triggering the im- 
mune response, which may result in variations in speed, 
timing, and magnitude of molecular responses to the 
pathogen [48,51]. 

Do patterns in transcriptomics distinguish Rps genes 
from alleles? 

An allele is defined as an alternative form of the same 
gene that is located at a specific position on a specific 
chromosome. In general, different alleles of the same 
gene are determined by genetic tests for allelism. How- 
ever, due to the complexity of phenotypes, such as vari- 
ation in the proportion of surviving plants of a pure NIL 
that contains an Rps gene, as observed in this study, 
such a test may not be effective in determining whether 
the resistance carried by two individual lines is controlled 
by different alleles or different genes. Moreover, many 
NBS-LRR gene models predicted in the Williams82 gen- 
ome are clustered in a small number of genomic regions 
[32,52], which makes it more difficult to distinguish genes 
resistant to the same pathogen by classic genetic analysis, 
especially if they are physically adjacent or closely linked 
to each other. 

The Rps genes/alleles in the 10 NILs have been previ- 
ously mapped to three genomic regions: Rpsl-a, Rpsl-b, 
Rpsl-c, and Rpsl-k on chromosome 3 [6,53,54], Rps3-a, 
Rps3-b and Rps3-c on chromosome 13 [1,54,55] and 
Rps4, RpsS and Rps6 on chromosome 18 [1,54,55]. How- 
ever, because limited numbers of molecular markers 
were available and used when these genes/alleles were 
identified, the genetic distances between genes or alleles 
mapped in the same chromosomal regions have not 
been determined. Actually, the alleles at the Rpsl and 
Rps3 loci were simply designated based on their disease 
reaction to a set of P. sojae races. As these three regions 
are highly enriched with NBS-LRR genes, according to 
the soybean reference genome, it is possible that some 
designated alleles of the same gene may be located at 
different loci. 
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As different alleles of the same gene are generally 
more identical to each other than to other genes, people 
are inclined to speculate that the defense response medi- 
ated by different alleles of a same gene may be more 
similar than those mediated by different genes. At first 
glance, this hypothesis seems to be supported by the ob- 
servation that the NILs with Rps3-a and Rps3-b exhib- 
ited similar patterns of transcriptional changes upon 
inoculation with P. sojae, as did the NILs with Rpsl-b 
and Rpsl-c (Figure 2). However, the pattern of transcrip- 
tional changes mediated by Rps3-c was found to be most 
similar to those mediated by Rps4, Rps5, and Rps6. Thus, 
the patterns of transcriptomic changes do not appear to 
predict genetic similarity and difference among these 
Rps genes and alleles, although we could not rule out 
the possibility that Rps3-c and Rpsl-a are non-allelic to 
the Rps3 and Rpsl loci, respectively. 



apparent that some of the genes belonging to these families 
play pivotal roles in defense signaling transduction. During 
compatible interaction, the avirulence effector of P. sojae 
induces a set of downstream responses that lead to disease 
development (Figure 6B). The signatures of the susceptible 
responses would include up-regulation of the JA pathway, 
suppression of ET pathway and with no significant changes 
in SA and BR pathways. It seems possible that the defense 
responses might be delayed in the compatible interaction 
due to insufficient activation of the expression of genes re- 
lated to ROS and phytoalexin biosynthesis, as well as 
MAPK signaling. Together, these data gained from a unique 
set of soybean NILs provide a comprehensive picture re- 
garding the molecular mechanisms underlying incompat- 
ible and compatible interaction between soybean and P. 
sojae, which shed light on the nature and timing of molecu- 
lar responses mediated by individual Rps genes. 



An integrated picture of molecular responses to P. sojae 

Based on analysis of putative defense- related genes identi- 
fied in this study and the patterns of their expressions 
upon inoculation with the pathogen, we propose a hypoth- 
esis that four major phytohormone signaling pathways are 
involved in defense responses (Figure 6A). During incom- 
patible interaction, the Rps protein recognizes the aviru- 
lence effector of P. sojae, initiating signaling transduction 
that involves the SA, JA, ET and BR pathways. The SA, ET, 
and BR pathways are activated, whereas the JA pathway is 
suppressed. Although the specific functions of individual 
genes encoding TFs belonging to the WRKY family, as well 
as those genes for reactive oxygen species (ROS) and MAP 
kinase (MAPK) cascades remain to be elucidated, it is 



Conclusions 

Comparative transcriptomic analysis of ten resistant NILs 
and susceptible recurrent line Williams with and without 
inoculation allowed us to identify DEGs associated with 
defense responses to P. sojae that were unique in each NIL 
and common among all these NILs, and thus depicted 
functional "fingerprints" of individual T^s-mediated resist- 
ance. Further analysis revealed multiple branches of 
putative ET, JA, ROS, and MAPK regulatory networks 
underlying the defense responses. Such responses exhibited 
dramatic variations among the soybean NILs containing 
distinct Rps genes/alleles. We propose different timing of 
individual i^s-mediated defense signaling to the same 
pathogen may be largely responsible for such variations. 
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Figure 6 An integrated picture for gene network of interactions between soybean and racel of Phytophthora sojae. A: Incompatible 
interaction. B: Compatible interaction. Red = up-regulation of a pathway. Purple = down-regulation of a pathway. Grey = no significant change in 
pathway. Yellow = up-regulation of genes. Blue = down-regulation of genes. Dotted lines indicate interactions of pathways inferred from 
Arabidopsis studies. 
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Methods 

Plant materials and isolates of R sojae 

The susceptible cultivar Williams and its ten NILs, each 
containing a unique Rps gene/allele, were used this study 
(Additional file 1). Each NIL was generated by back- 
crossing the donor for at least five generations with 
Williams. For inoculation tests, the isolate pmg(l)-3 
(pathotype racel) of P. sojae was used. 

Treatment of soybean materials 

Seeds of each soybean line were planted in sterilized 
sand in 10.1 cm clay pots and placed in a greenhouse. 
Approximately seven days after planting, at the VC 
growth stage (unifoliate leave fully expanded), the seed- 
lings for each line were separated into four groups, each 
containing -20 seedlings. Two groups were challenged 
with P. sojae using a standard hypototyl inoculation 
method [56], while the other two groups were mock in- 
oculated (wounded) without the pathogen. At 24 hpi, 
stems were harvested from one set of inoculated and 
one set of wounded seedlings by excising 2 to 3 cm 
across the wounded site and storing immediately in li- 
quid nitrogen. The seedlings in the remaining two 
groups were kept for evaluating symptom development, 
which was recorded 7 dpi. Experiments were performed 
a total of three times. 

RNA sequencing and analysis 

Total RNAs were isolated using the RNeasy Plant Mini 
Kit (Qiagen Inc., Valencia, CA) according to the manu- 
facturer s instructions in conjunction with DNase. The 
quality of total RNA was determined using a Nanodrop 
spectrometer (Thermo Fisher Scientific Inc., Wilming- 
ton, DE) and 1% formaldehyde gel electrophoresis. Total 
RNA samples were then sent to the Genomics Center at 
Purdue University for sequencing and 101 bp paired-end 
reads were generated with the Illumina HiScanSQ sys- 
tem. Since samples were a mixture of RNA from both 
soybean and P. sojae, to exclude the effect of P. sojae se- 
quences, all reads that aligned to the P. sojae genome se- 
quence were eliminated [57]. All remaining reads were 
then aligned to G. max reference genome (v8.0, http:// 
www.phytozome.net) using TopHat software [58] with 
parameters set to allow only one mismatch, and 30 and 
100,000 bp of the minimum and maximum intron length 
based on the current gene annotation. Moreover, only the 
uniquely mapped reads or fragments were kept for further 
analyses. The raw count of each gene was calculated by 
HTSeq (http://www-huber.embl.de/users/anders/HTSeq) 
with the "intersection_nonempty" mode, and preceded to 
edgeR package [59] in the R-Bioconductor tools. 

To detect variability, we estimated the over dispersion 
from 40 highly confident soybean housekeeping genes 
collected from previously reported papers (Additional 



file 10) [24,33,60-63], considering different soybean lines 
as replicates of control and inoculation groups. The data 
were then modeled using a Negative Binomial model in 
edgeR to identify differentially expressed genes after in- 
oculation using the dispersion estimated from the house- 
keeping genes as a common dispersion. Any genes with 
coverage read count less than one count per million of 
two paired samples were removed in later analyses. Dif- 
ferential expression between the inoculation group and 
the mock-inoculation group for each line was tested for 
false discovery rate (FDR) that was controlled at 0.05 
[64] using edgeR package [59]. 

Hierarchical cluster and heatmap analysis 

Pvclust and hclust were performed using the log 2 FC of 
the 5,806 IIGs with default parameters using distance 
method "correlation" for complete linkage clustering 
analysis. Pvclust provides two types of ^-values to assess 
the uncertainty for each cluster: approximately unbiased 
(AU) p-value and bootstrap probability (BP) value, via 
multi-scale bootstrap re-sampling [35]. The heatmap 
representing the expression intensity and direction was 
drawn using pheatmap R package with the distance 
method "euclidean" for both rows and columns [65]. 

Gene Ontology, heatmap and homologous gene 
functional analysis 

Annotations of GO terms were obtained from the AgriGO 
website based on the G. max model [66], and GO bio- 
logical process categories were used. These terms were 
manually classified into six broad functional groups, re- 
sponse to biotic or abiotic stress! 'biological regulation; 
'growth and development; 'transport; metabolism; and 
miscellaneous; while genes without GO annotations were 
grouped as unclassified'. For genes with multiple GO cat- 
egories, only one category was selected based on priority. 
Homologous genes were searched against annotated Ara- 
bidopsis gene protein database (The Arabidopsis Infor- 
mation Resouces, www.arabidopsis.org) using BLASTP to 
verify gene functions manually. 

Quantitative real-time PCR analysis 

The qRT-PCR was carried out using StepOnePlus™ Real 
Time PCR system (Applied Biosystems). The RNA sam- 
ples used as templates for qRT-PCR were the same as 
those used for RNA-Seq. The cons4 gene [60] was used 
as internal control for normalization of qRT-PCR data. 
For each gene, three experimental replicates were per- 
formed. Pearson correlations were calculated between 
RNA-Seq and qRT-PCR methods using six out of seven 
randomly selected genes for Log2 fold change of inocu- 
lated and mock-inoculated treatments in each soybean 
line (Additional file 3). 
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The transcriptome sequences presented in this study 
have been deposited in NCBI Gene Expression Omnibus 
under accession numbers GSE48524. 
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